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Abstract. 

We consider a semiclassical approximation, first derived by Heller and coworkers, 
for the time evolution of an originally gaussian wave packet in terms of complex 
trajectories. We also derive additional approximations replacing the complex 
trajectories by real ones. These yield three different semiclassical formulae involving 
different real trajectories. One of these formulae is Heller's thawed gaussian 
approximation. The other approximations are non-gaussian and may involve several 
trajectories determined by mixed initial-final conditions. These different formulae are 
tested for the cases of scattering by a hard wall, scattering by an attractive gaussian 
potential, and bound motion in a quartic oscillator. The formula with complex 
trajectories gives good results in all cases. The non-gaussian approximations with 
real trajectories work well in some cases, whereas the thawed gaussian works only in 
very simple situations. 



1. Introduction 

The Feynman propagator {xf\K{T)\xi) = {xf\exp {—iHT/h)\xi) can be interpreted 
as tlie time evolution of a wave function ip[xf,T) that is initially localized at Xj, 
ip{x,Q) = 6{x — Xi). In the semiclassical limit {xf\K{T)\xi) can be obtained from 
the classical trajectories connecting the initial point Xi to the final point Xf in time T. 
The formula which does this is known as the Van Vleck approximation. 

One is more likely, however, to be interested in the propagation of a smooth 
wavepacket, say ipo{x)j rather than an eigenfunction of the position or momentum 
operators. The straightforward way to accomplish this propagation, given the knowledge 
of K{T), is with an extra integration, thus 



ip{xf,T)= {xf\K{T)\xi) dxi{xi\ipo) 



(1.1) 
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In the semidassical limit this formula becomes 



^sc{xf,T)= {xf\K{T)\x.) 



VanVlcck dXi(Xi|V'o) • 



(1.2) 



One problem with this approach is that there are usually several classical 
trajectories going from Xi to xj and they are hard to find, the so-called root-search 
difficulty. We shall return to a discussion of this point in section 7. In the present 
paper, we shall accept the root-search problem and start from equation ()1.2|) . 

The main purpose of this paper is to derive various approximations for Eq. ()1.2|1 and 
to test them in simple problems. Hence we restrict ourselves here to systems with one 
degree of freedom. In principle, all the expressions that we obtain are easily generalized 
to multi-dimensional systems. In practice, of course, the difficulty of the numerical 
calculations increases very fast with the dimensionality, and it also depends on the 
actual approximation used. 

We shall evaluate the integral by the method of steepest descent which is also often 
called, somewhat inaccurately, stationary phase. We shall see that the stationary point 
is generally complex, giving rise to complex classical dynamics. Such a calculation was 
performed first by Heller and collaborators Pel87[ IHel88j . The result that they derived 
can be shown to be identical to our Eq. ()2.14p . They used it to calculate the motion 
of a wavepacket in the Morse potential. Their work does not seem to have received 
from physicists all the recognition that it deserves. Perhaps this was because there were 
two rather lengthy papers which contained many somewhat unfamiliar notions. Their 
first paper |Hel87j reached the result through heuristic arguments and the second one 
|Hel88j . though more direct, made many references to the first. In our paper, on the 
other hand, we shall present in section 2 a very simple and very short derivation of 
Heller's formula, Eq.(j2IIll). 

This basic complex approximation will be the starting point for all other 
developments of this paper. We shall show that it can be handled for relatively simple 
problems, but we shall also derive three subsequent approximations involving only real 
trajectories, which are different for each of the three cases. One case yields the well- 
known Heller Thawed Gaussian Approximation (TGA) |Hel75j . The other two give 
different results and we shall show that they can be quite accurate in some situations. We 
shall illustrate the application of the various semidassical formulae with two very simple 
examples (section 4) and two not so simple ones (sections 5 and 6). We hope that our 
derivation of the basic formula and our examples comparing the several approximations 
will stimulate new interest in semidassical formulae derived from complex trajectories. 

In order to simphfy our calculations we take for initial state |^o) a gaussian 
wavefunction with average position and momentum q and p respectively, and position 
uncertainty Ag = h/ \/2. This state is also the coherent state \z) of a harmonic oscillator 
of mass /i and frequency u. It is defined by 




e"" |0) 



(1.3) 
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with |0) the harmonic oscillator ground state and 



1 -P 



1 (q , ■p\ 



:i.4) 



In the above q, p, and are operators; q and p are real numbers; z is complex. The 
parameters 

b={h/^uj)^ and c = (hfiu)^ (1.5) 

define the length and momentum scales, respectively, and their product is h. The final 
form of Eq. ()1.2|l . our semidassical approximation for wavepacket propagation, is then 



{xf\K{T)\z)^ {xf\K{T)\x,) 



'j/ VanVleck 



(lxi{xi\z) = ip^c{xf,z]T) 



1.6) 



2. Approximation with complex trajectories 

Before we perform the integral in equation (jl.fij) by stationary phase, we recall the 
relations between the elements of the tangent matrix m and the second derivatives of the 
action function S{xf, xf, T) computed along the real classical trajectory connecting Xi to 
X/ in time T. Given a classical trajectory X(t), P{t) with X{Q) = Xi and X{T) = Xf, its 
tangent matrix m connects, in the linearized approximation, a small initial displacement 
6xi,6pi about the trajectory at t = to the propagated displacements Sxf,6pf at time 
T. The relation between m and the second derivatives of the action is 



b 



( 



5,: 



Sif 



C Sii 



c 1 
b S,f 

Sif 



\ 



J 



( ^\ 



/ 



11 



IP 



\mpq nipp J 



[2.1) 



where Sa = d'^S/dxi'^, Sif = Sfi = d'^S/dxidxj and Sff = d'^S/dxf'^. A derivation of 
this formula can be found, for instance, in sec. (2.6) of [BarOlj . Notice that we define m 
using the coherent state scales b and c. Inverting this equation we obtain 



c m. 



11 



b m, 



if 



IP 



b m, 



S 



c m. 



ff 



pp 



IP 



b m, 



[2.2) 



IP 



We shall write our final results in terms of elements of the tangent matrix m, rather 
than in terms of derivatives of the action. In this notation the Van Vleck propagator is 



{xf\K{T)\x,) 



1 



VanVleck 



JS/h-in/4 



by/2 



rrm, 



[2.3) 



IP 



For short times rugp is positive and the square root is well defined. For longer times 



m 



IP 



may become negative by going through zero. This happens at the focal points. 
At these points the Van Vleck formula diverges and must be replaced by a higher 
order approximation involving Airy functions. Sufficiently away from the foci the 
approximation becomes good again, as long as one replaces nigp by its modulus and 
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subtracts a phase 7i/2 for every focus encountered along the trajectory. We shall not 
write these so-called Morse phases explicitly because they will cancel out in our final 
formula. 

The overlap {xi\z) is given by 

(xi -qf\ f i 



{Xi\z) =n *b 2 exp I — j exp y-^p{xi - q/2) j . (2.4) 

In the semiclassical limit the propagated wavepacket ()1.6p is then 



with 

$(x;, X., T)='-[S + p{x, - q/2)] - ^^^^ - ztt/A . (2.6) 

We shall now approximate the integration over Xi by the stationary phase method. 
To do so we assume that the pre-factor in ()2.5p is a slowly- varying function of Xi, which 
means that the stationary point Xq will be computed by imposing zero variation on $ 
alone. The pre-factor will be simply computed at Xq and will not be expanded. We 
refer to section 3.3 of [BarOlj for a discussion of this procedure. We need to calculate 



the first and second derivatives of $. The first derivative is 



where we have defined 



dS 

Pi{xf,Xi;T) = -— . (2.8) 



The second derivative is 



9^$ i dvi 1 
$" = ilJL = _ _ . (2.9) 

dxi^ h dxi IP 

The stationary position xo{xf,T), solution of the stationary phase condition $' = 0, is 
given by the equation 

^ + z^ = (2.10) 

c 

where we have used c = h/b and we have made the definition 

Po{xf ; T) = pi{xf, Xq]T) . (2.11) 



Equation ()2.10|) makes it clear that Xq is usually complex, and therefore the stationary 
classical trajectory itself is complex. It leaves xq at time with momentum po, also 
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complex, and it arrives at Xf, which is real, at time T. It goes without saying that, for 
such a complex trajectory, the tangent matrix m is complex as well. 

To integrate over Xi, we expand $ around xq to second order, which means that 
we write $ ~ ^{xq) + $"(^o) i^i ~ ^o)^ /2 • If lm{dpi/dxi) is less than h/b"^ at Xq, then 
Re($) has a negative definite quadratic term and the integral can be performed in the 
gaussian approximation. Assuming this to be the case we find 



Finally, using equations (j2.9j) . (j2.8j) . (j2.2j) . we notice that 

$"(xo) = Is,. - ^ = l^!!hl±^ . (2.13) 
n 0^ 0^ mqp 

This simplifies the pre-factor in ()2.12|) and we get 

iJCT{Xf,Z]T) 



X exp 



^ S'(x/,xo ;T) + ^]9(xo - g/2) 



" ^' ' 262 



(2.14) 



We have changed the subscript on ip from "sc" to "CT" to indicate that this semiclassical 
approximation is calculated with complex trajectories. 

Equation ()2.14|1 . the semiclassical limit of the propagated wavepacket, is the Heller 
result |Hel87| IHel88j . It is the basic result from which all other approximations in this 
paper are derived. It is not an initial value representation and it may involve more 
than one complex classical trajectory. At first glance one might think that this formula 
represents a frozen gaussian, since the quadratic term in the exponent has a fixed width 
b. A closer look, however, reveals that, since the classical trajectory is complex, all other 
terms in the exponent also contribute a real part, which changes the effective width: 
the complex character of the stationary trajectory makes the wavepacket thaw and, at 
the same time, assume non-gaussian shapes. Notice that no Morse index or phases of 
7r/4 appear in ()2.14|1 . At T = one has rugq = 1, m^p = : one chooses the positive 
square root and the overlap ()2.4j) is recovered. For other times one simply follows the 
phase of the complex number rugg + inigp to get the phase of the pre-factor. 

Before we close this section we re-write the boundary conditions for the complex 
trajectories in a more convenient form. The stationary trajectory starts at xq (usually 
complex) and ends at x/ (always real). The initial momentum (also complex) is po, 
given by ()2.1H) and ()2.8|) . Since the coordinates X(t) and P(t) along the trajectory are 
both going to be complex, it is convenient to define the variables 

1 fx P\ 1 fx P\ , , 

where u v* in general. In this notation Hamilton's equations read 

dH , ^ dH , , 

inii = and — inv = —— . (2-16) 

ov ou 
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Finally the stationary condition (|2.1(jp can be re-written as 

— + i— = - + i- 2.17 

c b c 

which is equivalent to u{0) = z. Therefore the stationary trajectories are solutions of 
Hamilton's equations satisfying 

u{0) = z and X{T) = xj . (2.18) 



3. Approximations with real trajectories 

The stationary phase approximation of the previous section replaces the integral over 
a continuum of classical real trajectories by a few complex ones. Finding complex 
trajectories, however, is generally harder than finding real trajectories and one is 
tempted to look for further approximations to Eq. ()2.14j) in terms of real trajectories 
only. These approximations should be good if the stationary complex trajectory is 
sufficiently close to a real trajectory. 

What real trajectory? According to ()2.18j) . the complex trajectory is determined by 
two pieces of data: a final position Xf which is real and an initial coordinate u{0) = z 
which is complex; see fl2.17|l . The latter is made up of a real position q and a real 
momentum p. Thus we have a total of three real parameters, Xf, q , and p , that could 
be used to determine a real trajectory. Since it actually takes only two parameters, we 
have three obvious ways of choosing a real trajectory that might sometimes be a good 
approximation to the complex one. The first way is to give the two initial conditions 
{q,p), which are the center of the initial wavepacket in phase space. The other two 
ways use the final position Xf with a single one of the initial conditions, either {xf,q) 
or {xf,p). Each of these three choices leads to a possible approximation in terms of a 
real trajectory and, as we shall see, they are all different. In this section we shall derive 
formulae for each of these three cases. The extent of their validity will be examined in 
later sections. 

We can try to get a rough sense of how close these three real trajectories are to 
the unique complex one. The latter goes through Xf : hence the {xf,q) and {xf,p) 
trajectories have at least the merit of being at the right position for the final time T. 
To judge the closeness at time 0, we can look at formula ()2.14|) which tells us that, if 
Im(a;o) is sufficiently small compared to b, the semiclassical wave function at time T 
will be very small unless Re(xo) — g is of order b or smaller. This is a small distance 
in the semiclassical limit since b is of order h}^"^. We can also estimate the closeness of 
the initial momenta by looking at the stationary phase requirement ()2.im) . separating 
its real and imaginary parts thus 

Rexo — q Impo 

7"^ ~ ~r (3-1) 

tie po — p Lm Xq 
c b 
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Again assuming Im(xo) sufficiently small compared to b, we see from the second equation 
that Re(po) —p is of order c or smaller, and c is again proportional to h^^'^. The conclusion 
is that a semiclassical approximation in terms of real classical trajectories might be useful 
if conditions are right. 



3.1. The central trajectory approximation: Heller's TGA 

To begin, we choose the real trajectory defined by the two initial conditions {q,p), the 
central trajectory of the packet. We call {qx, Pt) the real end point of this trajectory at 
time T. This trajectory will now be the 0-order approximation. All values of S and its 
derivatives that occur in the equations, including the m matrix, will be taken for this 
trajectory, unless we state otherwise explicitly. All these quantities are now real. 
We write 

Xq = q + Axj (3.2) 

Xf = qr + Ax/ (3.3) 
dS 

Po = ~^(^/' xo;T) =p + Api=p- SiiAxi - SifAxf (3.4) 
dS 

Pf = +^(^/' xo;T) = pt + Apf =pt + SfiAxi + SffAxf . (3.5) 

The stationary phase condition ()2.1()|1 . can be re-written as Axi/b + iApi/c = 0. Then 
equation ()3.4|) can be solved for Axi 

= - o /. = ^^f (3-6) 

Sii + ic/b nigq + iniqp 

To simplify the notation, the complex number m^^ + inigp will be called m+ . 

We proceed to expand the exponent in ()2.14j) about the real trajectory. As before, 
the pre-factor is assumed to be a slowly-varying function and will be simply computed 
at the real trajectory. We write 

^S{xf,xo;T) + ^p{xo - g/2) - ^^^^yf^ 

■ -S - -pAxi + -prAxf 
n n n 

+ ^{SiiAx,^ + 2S^fAxiAxf + SffAxf^) 

+ ^M/2 + ^pAx.-^. (3.7) 



The linear terms in Axj cancel. The quadratic terms can all be written in terms of 
Axf = {xf — qr) using ()3.6|) . All second derivatives of the action can be written in 
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i 

— {SuAx,' + 2S,fAxAxf + SffAx/) - 



Axl 
" 262 



m 



qq 



2i- 



vv 



qp 



771: 



qp 



Axl 
262 

1 Axf' 
m+'^mqp 262 



[imqq - 2i{77lqq + im^p) + i77lpp7n+'^ - 77lqp] 

[~i7nqq + mqp + impp7n+'^] 



i Ax/ 

[1 - m+TRpp] 



Axf"^ mpp — im. 



pq 



m^mqp 262 



262 _|_ 



"qp 



where in the last equality we have used rnqqiripp — rnqpiripq 



With these simplifications the semiclassical wavepacket formula becomes 

5-1/2^-1/4 



7pqp{xf,z;T) 



niqq + i7n 



exp 



qp 



771. 



VP 



7771. 



pq 



qg 



7771 



qp 



Xf - Qt 



(3i 



(3.9) 



The subscript qp on indicates the variables used to compute the real trajectory. The 
spreading of the propagated wavepacket is now explicit in the coefficient of the gaussian 
term. This is exactly Heller's thawed gaussian approximation or TGA |Hel75j . As 
discussed at the end of section 7, a similar approximation [BarOlj can be obtained with 
coherent state path integrals. 



3.2. The approximation by trajectory q 



X 



f 



We consider now as 0-order the real trajectory which starts at q and ends a.t xj after 
time T. We call pi its initial momentum, which is a function of Xf, q, and T, and we 
write 



Xq = q + Ax. 
dS 



Po 



dxi 



(x/, xq]T) =pi + Api =pi- SiiAxi 



(3.10) 



Notice that the complete expansion of Pq to first order should be Pi — SnAxi — SijAxf 
but, as X/ is fixed, Axf = 0. Eq. ()2.10|) gives us the relation between Axi and Api 

■ (xo - q) 



Po-P _ Po-Pi _^ Pi-P 



b c c 

Thanks to Eqs.ljSHHD and ^TTT^ we obtain 



7C 



Api = —Axi - {pi -p) = -SiiAxi 
6 



(3.11) 



(3.12) 
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which gives, in terms of the tangent matrix, 

q .^c c mgg + irrigp 



We expand the exponent of Eg . (12.141) around this real trajectory. Once again we take 
the action S, its derivatives Si, Sj, and its second derivatives Sa, Sif, Sjf, including the 
m matrix, for the 0-order trajectory, unless we state otherwise explicitly. We obtain 

^:S{xf,Xo]T) + ^p{xo ~ q/2) - 



b + OjAxj + -DjjAxj + —pAxi + —zpq 



h V 2 " ' ; ' r ' ' 2/1"" 262 

+ 2^^^ + - ^^^^^^ + 2 I X - 6^ J 

i i irugp {pi -pY ^ inigp {pi - pf 
o + —rpq — ^'^ , L If 



(3.14) 



h 2h nigg + irUgp (? 2 jTj,^^ _|_ ivrigp (? 

Finally, the semiclassical propagator as a function of g, p, x/, and pj, becomes 



l5(x„g;T) + ^M 



2 + zmgp 



. (3.15) 



The gaussian in the exponent is now in the difference between pi, the initial momentum 
of the real trajectory, and p, the initial momentum of the center of the wave packet. 
Notice that there might be more than one trajectory satisfying the boundary conditions 
X(0) = g, X{T) = Xf . Notice also that the wave function is not restricted to a gaussian 
anymore. And that this is not an initial value formula. 

3.3. The approximation by trajectory p — > Xf 

In the same way, if we choose as 0-order the trajectory specified by p and Xf instead of 
q and Xf, we call the initial position for this trajectory, which is a function of Xf, p, 
and T, and we write 

Xo= qi + Axi 

dS (3.16) 

Po = ~^(^/' xo]T) =p + Api=p- SiiAxi . 

Thanks to eqs. ()2.10|) and ()3.16|) the new expression for Axj is 

Ax. = '^^^ (g - g.) . (3.17) 



Semiclassical Propagation of Wavepackets with Complex and Real Trajectories 



10 



In the end we obtain a propagating semiclassical wavepacket different from the other 
two, namely 



ipXfp{xf,z;T) 



5-1/2^-1/4 



exp 



qp 



ls{xf,qr,T) + ^^pq 



m, 



h I rUqq + iniqp 



The gaussian in the exponent is in the difference between gj, the initial position of the 
real trajectory, and q, the initial position of the center of the wavepacket. This formula 
shares many of the features of the previous one, including not being a gaussian and not 
being an initial value formula. 

The three real trajectory approximations we derived here are formally similar. 
Among the three real variables q, p, and xj, two are chosen to fix the trajectory. 
Then the corresponding formula carries a gaussian damping factor in the third variable. 
In the next three sections we discuss some examples that will help us compare these 
approximations with ipcT which involves complex trajectories. As usual all these 
formulae are exact for the free particle and the harmonic oscillator. 



4. First applications: the free particle and the hard wall 

The exact result for the propagation of a free particle wavepacket is of course well- 
known, but it can serve as a test of the formulae of sees. 2 and 3. And it is instructive 
to be able to see explicitly the different trajectories involved in each approximation. 
The wave function for the packet at time is given by ()2.4j) . The hamiltonian is simply 
H = where /i is the mass. The action is 

S{xj^xc,T) = ^^^^l^^. (4.1) 

All trajectories have constant momentum or velocity. In addition to the other 
parameters we shall use lo = cj = fij ^b"^, which is the frequency of the oscillator upon 
which the coherent states are built, and v = p/fi, the central velocity of the wavepacket. 
The elements of the tangent matrix are niqq = 1, niqp = ujT, ntpq = 0, nipp = 1. 

The exact expression for the propagated wavepacket follows from elementary 
quantum mechanics and can be written 



(xf — q — vT)'^ .if q vT' 



X exp 



262(1 + cj2T2 



' s i f q vT\ 

<^-'-^T) + -p^x,----j 



(4.2) 



It describes a gaussian packet of constant momentum p, whose real width byl + uj^T^ 
increases with time, but whose total width is actually complex and given by b\/l + iuT . 
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For the approximation of sec. 2, the stationary point Xq of the integration, which 
is the origin of the complex trajectory, and its associated momentum pq, which is the 
constant momentum of the trajectory, turn out to be 



Xf + iurL + i^p] p + i'l^x^^q) 

One can check the two boundary conditions ()2.18|) 

^ + ^— = 7 + «- xo H T = Xf . (4.4) 

b c b c fi 

One can also verify that, when Xf has the value q + vT, one finds the simple answers 
Xq = q and po = p : the trajectory is just the central trajectory of the packet. 

Let us look also at the three approximations with real trajectories. For the {q,p) 
case, formula ()3.9j) contains both and pt ■ The former is clearly qr = q + vT, while 
Pt = p, the constant momentum. For {xf,q), formula ()3.15|) . we need Pi which is 
Pi = fi{xf — q)/T. We don't need p/, but it is obviously the same as pi. Finally, for 
{xf,p), formula ()3.18|) . we need qi = xj — pT/fi, since the constant momentum is p. 
Once again, for the free particle, all four semiclassical approximations give the exact 
answer. 

To make things less trivial, we shall now have the wavepacket bounce elastically 
against a hard wall placed at the origin. The packet arrives from the right with a 
negative momentum and we assume that it starts far enough to have no appreciable 
value at the wall at time 0. The only relevant spatial region is x/ > ; the other side of 
the wall does not exist for this problem. Of course the wall could be an approximation 
to a very steep, but not totally hard, potential, in which case there would be some 
barrier penetration to the other side. And it would be interesting to see what happens 
in the limiting process, but we shall not do that here. 

The exact solution is quite simple, once one knows the exact free particle formula 
()4.2|) : you take the Xf < part of the latter, you reflect it with respect to the wall, and 
you give it an additional minus sign. Thus the complete solution has two parts, both 
restricted to x/ > 0, the original free particle part, and the reflected part which uses 
the X/ < piece of the free particle : 

V'wall , exact free, exact 

(x/,z;T) -^free,exact(-a;/,^;T) 

(4.5) 

X/ > only , T > 

This construction ensures that the Schroedinger equation is satisfied everywhere and 
that the total wave function vanishes at the wall. 

Next we look at the semiclassical approximation with complex trajectories, sec. 
2. Given Xj and x/, both > 0, there are two ways to go from one to the other in 
time T, the direct way and the way which bounces off the wall. The first trajectory is 
identical in all respects to a free particle's: it has the same action 1)4.11) and the process 
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of finding the stationary point xq of the integration over Xi is the same. Hence, for this 
direct trajectory, formula ()2.14|) yields the free particle result, restricted to Xf > . The 
reflected trajectory, on the other hand, has a different action. The distance between Xi 
and Xf, including reflection, is D = Xf + Xi, the speed is {xj + Xi)/T, the momentum 
fi{xf + Xi)/T, and the energy + Xj)^/2T^. Hence the action is 



Xf+Xi fi^Xf + XiY fl{Xf + 



\2 



S{xf, X,, T) = PD-ET = fi^-^{xj + X.) - I T = | . (4.6) 

It differs from ()4.1|) by a sign change for one of the two coordinates. The initial and 
final momenta follow from this action in the usual way: 

dS Xf + Xi dS 

Compared to the direct trajectory, the only change that needs to be made when we look 
for the stationary point of the Xi integration and we apply formula ()2.14|1 . is that the 
constant quantity Xf must be replaced everywhere by — x/. Thus we get again the free 
particle result, but for the value —Xf of the position. This approximation, therefore, 
yields the exact answer ()4.5|1 . except for the minus sign in front of the second term. This 
minus sign due to reflection is expected on very general grounds. It is a special case of 
the Morse phase mentioned in sec. 2, and it is understood easily if one looks at a soft 
barrier, such as a finite square step, and one goes to the infinite limit. 

Now let us try to apply to the hard wall the approximation with the real 
trajectory of sec. 3.1. We shall see that it does not work for this kind of problem. In 
this formulation there is a single trajectory, which begins at and evolves with T. 

Its endpoint is given hy qt = q — \v\T until it hits the wall, which happens for time 
= Q'/I^l- After this time qt becomes \v\T — q. Formula ()3.9|) yields for a given T a 
single gaussian wavepacket whose center follows this trajectory. The packet is incoming 
for T < Tj. and outgoing for T > T^.. Incoming and outgoing are never present at the 
same time and therefore there are no interference effects. In a calculation for a soft wall, 
one would see the full amplitude of the gaussian packet on the back side of the wall, 
with no damping due to barrier penetration, and again no interference effects. 

On the other hand, the other two approximations with real trajectories, (x/, q) and 
{xf,p), do work fine for the hard wall and give the exact answer. Each of them contains 
a direct and a reflected trajectory at all times, hence a sum of two wavepackets with 
interference. 



5. The inverted gaussian potential 

The hard wall is probably the simplest example after the free particle. However, except 
for ipqp, the calculation involves two trajectories. Our aim here is to consider a potential 
where a single classical trajectory contributes to the semiclassical formulae, to avoid 
dealing with interference and focus only on the differences arising from the complex or 
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real character of a trajectory. We choose an inverted gaussian potential with hamiltonian 

H{X,P)^^-e-''\ (5.1) 

The system has a single parameter, which is the width of the wavepacket. We set h ~ 1 
and we place the packet initially at g = —5, p = 1/2, hence central energy 1/8. There 
is a single bound level with energy -0.48, which excludes two simple limits: that in 
which the potential is just a perturbation, and that in which the problem is highly 
semiclassical. Here both the real and the complex trajectories need to be calculated 
numerically, as well as the exact packet propagation of course. The numerical methods 
for the trajectories will be discussed in the next section. 

Fig. 1 shows the wave packet at T = 7 for the different approximations, compared 
with the exact result, for two values of the width, b — 0.5 and b — 1.0 . The packet 
spreads very fast and rapidly becomes highly non-gaussian. It is remarkable how well 
this non-gaussian character is reproduced by some of the approximations, though not all. 
We shall have more to say in section 7 about comparing the various approximations. The 
top two graphs are the probability density versus Xf. The bottom two are the phase 
of the wave function. Both the modulus and the phase are very well reproduced by 
ipCT , the approximation via complex trajectories, and ipxfq , one of the real trajectories 
approximations, ip^fq being just as good as ipcT for this. On the other hand ipgp, which 
is a pure gaussian, does not do well at all. The small zig-zags in the approximate 
calculations of the phase are due to numerical imprecisions. Here, as well as in the 
next section, we do not show ip^fp , because it displays basically the same features as its 
counterpart ip^jq ■ Calculations with other values of b show these results to be robust, 
in the sense that ip^^g is always very similar to ipcT and that both agree well with the 
exact propagation. 

6. The quartic oscillator 

In this section we apply our semiclassical formulae to the case of a totally binding 
potential. We choose the quartic oscillator because it is probably the simplest system 
after the harmonic oscillator (for which all semiclassical formulae of sections 2 and 3 are 
exact). The hamiltonian is 

— + AX^ + BX^ (6.1) 

and the parameters are set to A = 0.5, B — 0.1, h — 1 . For these values the ground 
state energy is Eq a; 0.559 and the first two excited states have Ei ^ 1.770 and 
E2 ~ 3.319. For the wavepacket we choose q = 0, p = —2.0, and b = 1.0. This 
gives E = H{q,p) = 2.0 for the energy of the central trajectory, r ^ 4.7 for its period, 
and Xturn ~ ±1.6 for its turning points. Fig. 2 shows a sequence of four snapshots in 
the exact time evolution of the packet. The energy is low, the wavelength is large, and 
the interference effects are important. 
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As in the previous example, both the real and the complex trajectories have to 
be computed numerically. The real trajectories from q to xj can be calculated with a 
simple 'shooting method'. Since we need the propagated wave packet at all values of Xf, 
we simply integrate the equations of motion starting from X{0) = q (fixed) and several 
values of the initial momentum P(0) = pi. For each pi we obtain a final coordinate 
X{T) = Xf at which we evaluate the wave function. 

The calculation of the complex trajectories is more involved and we shall describe 
it in some detail. We follow the procedure introduced by Klauder |Kla87j (see 
also |Hel87| IAda89t IRub95j and |Rib04j for a different approach), which consists in 
propagating trajectories starting from 

X{0) = q + w P{0)=p + i'^w (6.2) 

where w = a + if3 is a complex number. The first of conditions ()2.18p is automatically 
satisfied for all w. The search method consists in propagating trajectories for all possible 
w's and picking those satisfying the second condition ()2.18p . X{T) = Xf . 

As discussed in |R,ub95j . fixing q, p and T and integrating Hamilton's equations 
with initial conditions ()6.2j) gives, for each w, a final coordinate X{T) = Xt, usually 
complex. This can be viewed as a map Xt = Xt{w). The values of w we need are given 
by the inverse map of the real X axis. If the map is analytic at w, then it is conformal. 
This implies, among other things, that the map is one-to-one in the neighborhood of w. 
If, on the other hand, w is a critical point of the map, a richer structure develops. If, for 
instance, dXf/dw = 0, but the second derivative is non-zero, it can be shown that the 
map becomes two-to-one in the vicinity of w, meaning that two different trajectories 
(corresponding to two distinct initial conditions) satisfy the same boundary conditions. 
This is the basic mechanism that generates multiple trajectories in systems with one 
degree of freedom. For short times the critical points of the map generically lie very far 
from the origin w = 0, corresponding to complex trajectories whose actions have large 
imaginary parts. As time increases, more and more of these points approach the origin, 
giving rise to significant contributions to the propagator |Rub95j . The singularities 
in the Xt{w) map produced by the critical points are called Phase Space Caustics 
(PSC). At these points the square- root in the pre-factor of Eq. ()2.14|l goes to zero and 
the semiclassical approximation fails (it actually fails in a finite neighborhood of the 
PSC's). 

For fixed g,p,andT, the complex trajectories form continuous families 
parameterized by Xf . Among the many families that might contribute to the 
semiclassical evolution of the wave-packet, one is special. This is the family that contains 
the real trajectory that starts at {q,p), and we call it the main family. If (qTyPr) is the 
end point of this real trajectory, then for x/ = qr the trajectory in the main family is 
real. In terms of the map, the point w = is mapped into Xj- = qj-. In some special 
situations, the main family alone may provide sufficient information for the semiclassical 
calculation jHel87| lHel02| IHel03] . but that is not always the case |Shu95t IShu96t lRib04j . 
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Fig. 3 (a) shows the map Xt = Xxiw) for T = 6.5. The hnes correspond to constant 
values of Re(X) and Im(X). The conformal property of the map guarantees that these 
hnes cross at right angle (provided the same scale is used for a and f3). The circle 
indicates the location of a critical point, corresponding to a singularity in the regular 
pattern of the lines. The thick lines correspond to families of trajectories for which 
Im(X(T)) = . The main family can be identified by its containing the point w = . A 
more careful look at this figure reveals the existence of several other singularities in the 
otherwise regular pattern of crossing lines. Each of them is a critical point of the map. 
The trajectories in their neighborhood, however, do not contribute significantly to the 
propagated wave packet, and we shall not take them into account. 

Let 

F = S{xf, xo; T) + pixo - q/2) + (6.3) 

be the exponent in Eq. ()2.14|) . For real trajectories Im(F) is > 0. Fig. 3(b) shows a 
greyscale topographic map of Im(F) for the trajectories calculated in Fig. 3 (a) and in 
the same (a,/3) plane. The continuous thick line is the main family seen on Fig.3(a). 
The line for the other family, which we call secondary, is part thick and part thin. The 
reason for this distinction is the following. The imaginary part of F for trajectories in 
the main family can be seen to be always positive. It starts with Im(F) = +oo for 
Xf = +00, decreases to zero at x/ = gr and grows to infinity again as Xf — oo. For 
the trajectories in the secondary family, however, lm{F) starts at +oo for x/ = +oo, 
and decreases steadily to Im(F) = — oo at x/ = — oo. This means that this family 
cannot be included in the semiclassical calculation for all values of x/. There must be 
a cutting point after which the family cannot be considered, otherwise it would cause 
the wave function to diverge. This phenomenon has become known as the problem of 
non- contributing trajectories, and has been studied by Adachi |Ada89j . Berry |Ber89j . 
Klauder |Rub95j and others |Shu95t IShu96t ITan98| IRib04j . The part of the secondary 
family shown with the thin line has not been included in our calculations. We note that, 
in |Hel87j . Huber and Heller do a similar calculation for a wave packet moving in the 
Morse potential. 

To understand qualitatively the role of secondary families in the semiclassical limit, 
we consider h to be really very small. In this case the (x/, q) and {xf,p) real trajectory 
approximations of section 3 become exact, and only the main family (actually only a 
small neighborhood of the real trajectory) contributes significantly. Demanding the 
approximation to be uniformly valid as ^ ^ 0, and using the subscript m for the main 
family and s for a secondary one, we must apply the following rules: (a) Trajectories with 
Im(Fs) < should be removed |R,ub95j . This avoids the divergence of exp {lm{Fs) / h} 
in the limit — > 0; (b) Trajectories with Im(Fs)(x/) > but lm{Fs{x j)) < Im(Fm(x/)) 
should also be removed to guarantee that, as ^ — 0, the main family always gives the 
dominant contribution; and (c) The discontinuity introduced by the sudden removal of 
a secondary contribution should be minimized. This criterion ultimately determines the 
choice of the cutoff point, which is located on the so-called Stokes line jR,ub95j . 
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Fig.4 shows the comparison between the exact calculation (thick line) and the 
different approximations of sections 2 and 3 for T = 2.5 and T = 4.5 . The most 
interesting feature here is that the real trajectory formula ipxfq becomes discontinuous. 
This can be understood with the help of panels (b) and (d) which show the initial 
momentum pi of the contributing trajectories as a function of x/ (for fixed q and T). 
We see that several branches appear, but it turns out that only the one shown with a 
thick line contributes significantly. Since this main branch covers only a finite range of 
Xf, the wave function drops to zero suddenly. In fact, the wave function diverges at the 
ends of the branches. In our calculations we have cut the branch a little before its end 
points. 

Fig.5(a) shows the real part of ipcri^f, z; T) for T = 6.5 calculated with the separate 
contribution of each of the families shown in Figs. 3. The abrupt cutoff of the secondary 
family is clear in this figure. Fig.5(b) shows again the separate contributions of the same 
families to \ipcT{xf,z;T)\'^ (dashed and solid lines) and their combined contribution 
(thick solid line). Notice the interference between the two families producing the 
oscillation in the probability density. 

Finally Fig. 6 shows \ipcT{xf, z; T) p, the calculation with complex trajectories. Each 
snapshot shows the exact result (thin solid line), the thawed gaussian approximation 
ipqp (dashed hne) and the complex trajectories approximation (thick line). We do not 
show ipxfq in these plots because the approximation is not good. The improvement 
obtained with the complex trajectories is clear. For T = 6.5 and T = 8.5, however, 
we can see a spurious peak close to the turning points. This might seem to result 
from cutting off the secondary family at the wrong point, as suggested in |Hel88j . A 
careful analysis shows that it is not the case: cutting off the secondary family at a 
smaller value of Xf would produce a sudden dip in the probability density followed by 
another sudden increase, as can be seen from Fig. 5(b). The spurious peaks are due 
to the phase space caustic that shows up close to the turning points. For the present 
situation, which is not quite in the semiclassical domain, the caustics have a strong effect 
on the semiclassical propagations and further corrections are necessary to improve the 
approximation. A uniform semiclassical approximation involving Airy functions can be 
derived by considering cubic terms in the saddle point approximation performed in the 
integral ()2.5j) . The derivation of this improved formula will be published elsewhere. 

7. Summary and discussion 

We begin by summarizing the results. We have derived and written down in this paper 
(sections 2 and 3) four semiclassical expressions for the propagation of an initially 
gaussian wavepacket, ^ipcT, i^qp , i^xfq , i'xfp ■ All four give exact results for the free particle 
and the harmonic oscillator. All except ipqp also give the exact result for the scattering 
by a hard wall. But ipqp is very bad for the hard wall: it produces no interference. We 
tested the formulae further on two examples of smooth potentials. For the attractive 
gaussian potential, the exact packet after some time is highly non-gaussian. Both ipcT 
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and i/jr^^g involve a single trajectory for each Xf, and both give very good results for the 
modulus as well as the phase of the wave packet. The same is true of ip^fp ■ But ipgp , 
which is a pure gaussian, gives a very bad fit. 

The other smooth potential is a quartic oscillator. Now there may be contributions 
from multiple families of trajectories and the results are quite different. We found in this 
case that the complex trajectories formula ipcT is superior to both the real trajectories 
approximation ^/'^^g and the single trajectory approximation ipgp. One reason is that 
the real trajectories contributing to ipxfq have finite branches, which causes a drastic 
discontinuity in the approximate wave function. The complex trajectories, on the other 
hand, form continuous branches that prolong into the complex domain. The projection 
of these complex families onto the real plane show that they follow closely the main 
branch of the real trajectories, but they continue to exist after the latter end. 

Thus the paper's main finding is that the complex trajectory approximation, ipcT 
of eq. ()2.14|) . gives very good results, at least for all the cases that we looked at. And 
one may well wonder how it could be that a relatively simple approximation, based on 
a single trajectory, works so well. The short answer is that, actually, the approximation 
is not that simple and it is not based on a single trajectory. But we shall go into more 
details. 

Many other approaches have been proposed for the propagation of an initial 
wave function, some of them especially concerned with the propagation in systems 
with many degrees of freedom. Among them. Initial Value Representations or IVR's 
|Mi1174l IHer84[ IHe191l |Kay94al IPolOai IThon4j have become especially popular . They 
are integral expressions over the phase space of initial conditions, in a spirit which tends 
to rejoin that of path integrals. The multiple integration is usually handled by Monte 
Carlo techniques |Kay94b| IZhaOSj . The initial impetus for these approaches jMill7nj 
was the root-search problem, i.e. the fact that the Van Vleck expression ()2.3|) requires 
one to find trajectories determined by mixed initial-final boundary conditions (starting 
at Xi and ending at Xf). This may lead to a difficult search in complicated situations, 
and in addition the root-search problem usually has several solutions, a fact which 
complicates further the time-development of the wave function. In the IVR's, on the 
other hand, the trajectories are determined solely by initial boundary conditions, as 
the name implies, and for each set of initial conditions there is only one trajectory, the 
unique solution of Hamilton's equations for these conditions. The work of getting the 
correct propagated wave function is done entirely by the integration, as is the case for 
path integral expressions. 

The present work is based on the Van Vleck expression and it is affected by root- 
search, except for the thawed gaussian formula ipgp of section 3.1. Root-search can 
truly be a burden, especially for multi-dimensional wave functions. More important, 
however, is the fact that characterizing the present approach as a one-trajectory or a 
few-trajectories approximation is highly misleading. We have used such characterization 
several times above, for instance in connection with the hard wall, or in the first few 
lines of section 5. But the truth is that there is a different classical trajectory for every 
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final Xf, and therefore tlie number of trajectories is really infinite. The work of getting 
the correct propagated wave function is done by them. 

The one case where we truly have a one-trajectory expression is ipqp ■ There, the 
same trajectory is used for all values of Xf. This case is truly an IVR, though a very 
simplified one, not involving any integration over the initial parameters. The only 
initial parameters used are those of the center of the packet. And it is also the case for 
which, in our explorations, the agreement with the exact result is worst, showing that a 
single classical trajectory is just not enough. Somehow, a variety of classical trajectories 
must enter into the formalism. In other words, we see the path integral idea reappear, 
Feynman's original view of the connection between quantum and classical mechanics. 

Among the various IVR's, there is one that has been especially successful and 
has had many applications, the Herman-Kluk approximation |Her84[ |Kay94a[ KTarOOl 
IWanOH INoiORj . It would be particularly interesting to have a direct application of 
the HK formula to the hard wall, since this is the simplest non-obvious case for which 
our complex trajectories formula is exact. We are not aware of the existence of such 
an application of HK, and it seems to us that the calculation will not be totally 
straightforward. 

Perhaps it is not superfluous to repeat once again that the main new result of 
this paper is that semiclassical approximations using complex trajectories, especially 
non-IVR approximations, are capable of being surprisingly good. All our semiclassical 
expressions generalize easily to many dimensions. The detailed presentation of these 
generalizations will be pubhshed separately, together with numerical applications. We 
note that two-dimensional semiclassical calculations with complex trajectories have 
already been considered in the context of coherent states. In |Hein2[ lRibn4j the quantity 
evaluated was the return probability amplitude, also known as the auto-correlation 
function. In |Shu95t IShu96| lUniOlj the authors, who seem to have been unaware of 
references |Hel87| lHel88j . calculate the propagated wave-function in the momentum or 
coordinate representation. But they do it for quantum maps, which is rather different 
from what we have done here since, in their case, there are no continuous classical 
trajectories, no Hamilton equations, and no Schroedinger equation. In the present paper 
we provided formulae for the calculation of the propagated wave-function for continuous 
systems in the coordinate representation. The basic equation for such calculations is 
the complex trajectory formula ()2.14j) first derived by Heller and coworkers in |Hel88j . 
We showed that one can approximate the complex trajectories by real ones, and that 
this can be done in various ways (we mention three of them), with different degrees 
of difficulty and of accuracy. The various formulae so generated will undoubtedly have 
different domains of applicability. And the complex formula, if it is practical, will always 
remain the better one. 

As our final topic, we shall present an alternative derivation of a formula very similar 
to ()2.14|) . which possesses the same three approximations in terms of real trajectories. 
The quantity we want to calculate is {xf\K{T)\z). To obtain the approximation ()l.(i|l . 
we introduced between K{T) and \z) the resolution of unity in terms of eigenstates of x, 



Semiclassical Propagation of Wavepackets with Complex and Real Trajectories 



19 



and then we replaced the exact x-propagator by its Van Vleck approximation. Instead 
of this, we shall introduce between {xf\ and K{T) the resolution of unity in terms of 
coherent states, thus 

{xj\KiT)\z) = J {xf\z')^{z'\KiT)\z) (7.1) 

and then we shall replace the exact coherent state propagator by its semiclassical 
approximation. The latter was worked out in detail in [BarOlj . Performing the integrals 
over q' and p' by the stationary exponent approximation, we arrive at an expression 
which is very similar, but not identical, to ()2.14j) . There are two important differences. 
One is that, in this case, the dynamics are governed by the smoothed hamiltonian 
T^ilyP) = {^\H\z) instead of the classical hamiltonian. The other is that the exponent 
acquires an extra term given by [Barnij 

This term plays a very important role and cannot be dropped from the approximation. 
Both 1 and the Ti dynamics get carried over into all the real trajectories approximations. 

Our main reason for mentioning this alternative is to stress that semiclassical ap- 
proximations are never unique and that there are always many ways to base an approx- 
imate quantal expression on a classical solution. Discussions of this appear in |Kla85j . 



who see it as a consequence of the over-completeness of the coherent-state basis, and 
also in IBarOlj . It is only in the limit of h going to zero that they all become the same. 
From our alternative to ()2.14|) involving the smoothed Hamiltonian and the correction 
term X, one can derive again three expressions in terms of real trajectories, different 
from those in section 3. Numerically, the two points of view yield roughly equally good 
(or sometimes equally bad) results. And once again, for quadratic Hamiltonians, all the 
approximations we have mentioned are exact. While in section 3 the ipqp approximation 
was the same as Heller's TGA, for the alternative formulation it is the same as the 
approximation in section 4 of |Bar01j . The latter reference contains a discussion of the 
differences between the two results. Both are initial value formulae using the real central 
trajectory of the packet, and both are gaussian in shape for all potentials and all times. 
Both are expected to fail if the potential is not very smooth on the scale of the size of 
the packet. 
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Figure 1. Parts (a) and (b) show the probability density in the inverted Gaussian 
potential for p = 0.5 and T = 7.0. In (a) b = 0.5 and in (b) 6 = 1.0. The thick solid 
line shows the exact result, the thin solid line correspond to ■0CT and the dotted line 
to ijjxjq (which cannot be distinguished from '0CT at this scale). Part (b) also shows a 
comparison with ipgp (symmetric Gaussian curve). Parts (c) and (d) show the phase 
of the wave packet, in units of tt, corresponding to (a) and (b) respectively. In this 
figure and those following, the x variable in the abscissa is actually the variable Xf of 
the text. 
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Figure 2. Exact quantum mechanical propagation in the quartic potential with 
A = 0.5, B = 0.1 and h = 1. The wavepacket is initially centered at g = 0, p = —2 
and has width b — I. The curves show the probability density at times (a) T = 
(thick line) and T = 2.5, (b) T = 4.5, (c) T = 6.5 and (d) T = 8.5. The period of the 
classical orbit of the center of the packet is r « 4.7. 
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Figure 3. (a) Map Xt ~ Xt(w) for T = 6.5. The lines correspond to constant 
values of Re(X) and Im(X) and the circle indicates the singularity. The thick lines 
correspond to the trajectories satisfying Iin(X7') = 0; (b) Gray scale topographic plot 
of the imaginary part of the exponent F for all the trajectories in (a). The shades of 
gray go from — oo (black) to +oo (white). The main family has Im(i^ ) > whereas 
the secondary family has a section where Im(F) < (thin line). 
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Figure 4. Parts (a) and (c) show the probabihty density in the quartic potential for 
T = 2.5 and T = 4.5 respectively. The thin solid line shows the exact result, the dotted 

line corresponds to ijjxfq, the dashed line to ijjqp and the thick solid line to ijjcT- Parts 
(b) and (d) show the initial momentum pi as a function of x for the {q, x) real trajectory 
for T = 2.5 and T = 4.5. The thick line shows the branch used in the calculation of 
ipxfQ and the star represents the central trajectory starting from q.p. The dashed line 
corresponds to the projection of the main family of complex trajectories into this real 
plane. 
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Figure 5. Separate contributions of the main (solid line) and secondary (dashed line) 
families for the wave packet at T = 6.5. (a) Real part of tp; (b)probability density 
Notice the abrupt cutoff of the secondary family. The thick line in (b) shows the total 
probability density, displaying the interference between the individual contributions. 




Figure 6. Probability density for the quartic potential for (a) T — 2.5; (b) T — 4.5; 
(c) T = 6.5 and (d) T = 8.5. The thin sohd line shows the exact result, the dashed 
line corresponds to Vgp and the thick solid line to Vct- 



